﻿using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Linq;
using System.Reflection;
using ESRI.ArcGIS.ConversionTools;
using ESRI.ArcGIS.DataManagementTools;
using ESRI.ArcGIS.DataSourcesGDB;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.GeoAnalyst;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geoprocessing;
using ESRI.ArcGIS.Geoprocessor;
using ESRI.ArcGIS.SpatialAnalyst;
using ESRI.ArcGIS.SpatialAnalystTools;

namespace Bes.Swsp.ToolBox
{
  /// <summary>
  /// Class for creating, loading, and saving StormDrainNetworks using ArcGIS data.
  /// </summary>
  public class StormDrainNetworkLoader
  {
    private StormDrainNetwork _stormNetwork;

    private IFeatureClass _nodesClass;
    private IFeatureClass _conveyanceClass;

    private IGeometricNetwork _geometricNetwork;
    private IFeatureClass _catchmentClass;
    private IRaster _flowDirRaster;
    private IRaster _runoffRaster;
    private IRaster _localAccumRaster;

    private bool _extractTerminalFlow = false;
    private string _flowAccumField = null;
    private string _loadAccumField = null;
    private string _tagField = null;

    // Storm Drain Network Field Names
    private const string ID_FIELD = "Id";
    private const string TYPE_FIELD = "Type";
    private const string DOWNNODE_FIELD = "DownNode";
    private const string UPNODE_FIELD = "UpNode";
    private const string UPNODES_FIELD = "UpNodes";
    private const string ENABLED_FIELD = "Enabled";
    private const string ERRORS_FIELD = "Errors";
    private const string SUBSURFACE_FIELD = "Subsurface";
    private const string SOURCECLASS_FIELD = "SourceClass";
    private const string SOURCEOID_FIELD = "SourceOID";
    private const string GLOBALID_FIELD = "GlobalID";
    private const string FLOW_ACC_FIELD = "accumulated_area_sq_ft";

    private Dictionary<string, string> _classNames = new Dictionary<string, string>();

    /// <summary>
    /// Gets or sets the storm drain network.
    /// </summary>
    public StormDrainNetwork StormNetwork
    {
      get { return _stormNetwork; }
      set { _stormNetwork = value; }
    }

    /// <summary>
    /// Gets or sets the feature class that represents the nodes in the StormDrainNetwork.
    /// </summary>
    public IFeatureClass StormNetworkNodeClass
    {
      get { return _nodesClass; }
      set { _nodesClass = value; }
    }

    /// <summary>
    /// Gets or sets the feature class that represents the conveyances in the StormDrainNetwork.
    /// </summary>
    public IFeatureClass StormNetworkConveyanceClass
    {
      get { return _conveyanceClass; }
      set { _conveyanceClass = value; }
    }

    /// <summary>
    /// Gets or sets the geometric network used to create the StormDrainNetwork.
    /// </summary>
    public IGeometricNetwork PipeNetwork
    {
      get { return _geometricNetwork; }
      set { _geometricNetwork = value; }
    }

    /// <summary>
    /// Gets or sets the catchment feature class corresponding to the local drainage area of each inlet in the pipe network.
    /// </summary>
    public IFeatureClass CatchmentClass
    {
      get { return _catchmentClass; }
      set { _catchmentClass = value; }
    }

    /// <summary>
    /// Gets or sets the flow direction raster.
    /// </summary>
    public IRaster FlowDirRaster
    {
      get { return _flowDirRaster; }
      set { _flowDirRaster = value; }
    }

    /// <summary>
    /// Gets or sets the runoff coefficient raster that represents runoff from each cell.
    /// </summary>
    public IRaster RunoffRaster
    {
      get { return _runoffRaster; }
      set { _runoffRaster = value; }
    }

    /// <summary>
    /// Gets or sets the Local Accumulation raster. This is the local surface accumulation without regard to flow through pipes.
    /// </summary>
    /// <remarks>This is normally calculated from the flow direction and runoff rasters, and does not need to be set. Only set if it has been
    /// precalculated prior to calling CreateStormNetwork(). If it is not set, it will automatically be calculated. Otherwise the set value
    /// will be used.</remarks>
    public IRaster LocalAccumulation
    {
      get { return _localAccumRaster; }
      set { _localAccumRaster = value; }
    }

    /// <summary>
    /// Gets or sets whether to extract overland flow lines for the last downstream outlets flowing to masked out rivers or to the edge of the grid.
    /// </summary>
    /// <remarks>Setting this value to true can significantly increase processing time, depending on input data (mainly the size of the flow direction grid).</remarks>
    public bool ExtractTerminalFlow
    {
      get { return _extractTerminalFlow; }
      set { _extractTerminalFlow = value; }
    }

    /// <summary>
    /// Gets or sets the name of the field that will store flow accumulation values
    /// </summary>
    public string FlowAccumulationField
    {
      get { return _flowAccumField; }
      set { _flowAccumField = value; }
    }

    public string LoadAccumulationField
    {
      get { return _loadAccumField; }
      set { _loadAccumField = value; }
    }

    /// <summary>
    /// Gets or sets the name of the field that will store tag values. Field must already exist.
    /// </summary>
    /// <remarks>Assumes tag values are correctly formatted (of the same type) as the given field.</remarks>
    public string TagField
    {
      get { return _tagField; }
      set { _tagField = value; }
    }

    /// <summary>
    /// Progress reporting event of operations performed by this class.
    /// </summary>
    public event EventHandler<ProgressEventArgs> Progress;

    /// <summary>
    /// Warning report for non-critical errors found during processing.
    /// </summary>
    public event EventHandler<WarningEventArgs> Warning;

    public void CreateStormNetwork(IWorkspace workspace, string name, int threshold, double flowConversionFactor = 1.0)
    {
      _LogStormNetworkStart(threshold, flowConversionFactor);

      _stormNetwork = new StormDrainNetwork();
      _InitializeStormNetwork(_stormNetwork);

      if(string.IsNullOrWhiteSpace(this.FlowAccumulationField))
        this.FlowAccumulationField = FLOW_ACC_FIELD;

      if (_localAccumRaster == null)
        _localAccumRaster = _CalculateLocalFlowAccumulation();
      else
        OnProgress(new ProgressEventArgs(1, 1, "Using existing local flow accumulation raster."));

      _ExtractLocalAccumulation(_stormNetwork.Catchments, _localAccumRaster);
      _stormNetwork.Accumulate();
      _stormNetwork.DetectErrors();
      _SaveStormNetwork(workspace, "overlapping_network");

      IFeatureClass streamClass = _CalculateTotalFlowAccumulation(_localAccumRaster, threshold);
      _RebuildSurfaceNetwork(streamClass);
      _stormNetwork.DetectErrors();
      _stormNetwork.DetectLowFlow(threshold);
      _SaveStormNetwork(workspace, name, flowConversionFactor);

      _LogStormNetworkEnd();
    }

    private void _LogStormNetworkStart(int threshold, double flowConversionFactor)
    {
      string pipeNetworkPath = GeoprocessingUtility.GetFeatureClassPath(this.PipeNetwork.OrphanJunctionFeatureClass);
      pipeNetworkPath = pipeNetworkPath.Substring(0, pipeNetworkPath.Length - 10);

      Logger.Log("---------------------------------------");
      Logger.Log(string.Format("Create Storm Network initiated by {0}", Environment.UserName));
      Logger.Log(string.Format("INPUTS:"));
      Logger.Log(string.Format("Pipe Network: {0}", pipeNetworkPath));
      Logger.Log(string.Format("Catchments: {0}", GeoprocessingUtility.GetFeatureClassPath(this.CatchmentClass)));
      Logger.Log(string.Format("Flow Direction Raster: {0}", GeoprocessingUtility.GetRasterPath(this.FlowDirRaster)));
      Logger.Log(string.Format("Runoff Coefficient Raster: {0}", GeoprocessingUtility.GetRasterPath(this.RunoffRaster)));
      Logger.Log(string.Format("Precalculated local accumulation: {0}", GeoprocessingUtility.GetRasterPath(this.LocalAccumulation)));
      Logger.Log(string.Format("Extracting Terminal Flow: {0}", this.ExtractTerminalFlow));
      Logger.Log(string.Format("Flow Threshold: {0}", threshold));
      Logger.Log(string.Format("Flow Conversion Factor: {0}", flowConversionFactor));
      Logger.Log("---------------------------------------");
    }

    private void _LogStormNetworkEnd()
    {
      Logger.Log("---------------------------------------");
      Logger.Log(string.Format("Create Storm Network completed."));
      Logger.Log(string.Format("OUTPUTS:"));
      Logger.Log(string.Format("Storm Network Nodes: {0}", GeoprocessingUtility.GetFeatureClassPath(this.StormNetworkNodeClass)));
      Logger.Log(string.Format("Storm Network Conveyances: {0}", GeoprocessingUtility.GetFeatureClassPath(this.StormNetworkConveyanceClass)));
      Logger.Log("---------------------------------------");
    }

    public void LoadStormNetwork()
    {
      _stormNetwork = new StormDrainNetwork();
      _LoadStormNetworkNodes(_stormNetwork);

      // We load catchments without shapes to allow for local flow accumulation calculations
      Dictionary<int, Catchment> catchments = _ReadCatchments(false);
      _LinkCatchmentsToInlets(catchments, _stormNetwork.Inlets);

      _LoadStormNetworkConveyance(_stormNetwork);
      _stormNetwork.DetectErrors();
    }

    public void LoadTotalFlowAccumulation(string flowAccumulationField, double conversionFactor = 1.0)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Reading accumulated flow values...", true), true);

      if (_nodesClass == null)
        throw new InvalidOperationException("The storm network node feature class is not specified.");

      if (_stormNetwork == null)
        throw new InvalidOperationException("The storm drain network is not loaded or created.");

      int flowField = _nodesClass.FindField(flowAccumulationField);
      if (flowField == -1)
        throw new ArgumentException("The specified flow field cannot be found in the storm network node feature class.");

      int idField = _nodesClass.FindField("Id");
      if (idField == -1)
        throw new InvalidOperationException("The specified storm network node feature class is invalid. It does not have an 'Id' field.");

      int progress = 0;
      int total = _nodesClass.FeatureCount(null);
      int increment = _CalcIncrement(total);

      IFeatureCursor cursor = null;
      try
      {
        cursor = _nodesClass.Search(null, false);
        IFeature node = cursor.NextFeature();
        while (node != null)
        {
          try
          {
            int id = (int)node.get_Value(idField);
            double flow = _SafeConvertToDouble(node.get_Value(flowField), 0.0);
            if (_stormNetwork.Nodes.Contains(id))
            {
              _stormNetwork.Nodes[id].TotalAccumulatedFlow = flow * conversionFactor;
            }
          }
          finally
          {
            ComUtility.ReleaseComObject(node);
          }

          if (++progress % increment == 0)
          {
            ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading accumulated flow values...", true);
            OnProgress(eventArgs);
            if (eventArgs.Cancel)
              throw new UserCanceledException();
          }

          node = cursor.NextFeature();
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} accumulated flow values read.", progress)), true);
    }

    public void LoadTotalLoadAccumulation(string loadAccumulationField, double conversionFactor = 1.0)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Reading accumulated load values...", true), true);

      if (_nodesClass == null)
        throw new InvalidOperationException("The storm network node feature class is not specified.");

      if (_stormNetwork == null)
        throw new InvalidOperationException("The storm drain network is not loaded or created.");

      int loadField = _nodesClass.FindField(loadAccumulationField);
      if (loadField == -1)
        throw new ArgumentException("The specified load field cannot be found in the storm network node feature class.");

      int idField = _nodesClass.FindField("Id");
      if (idField == -1)
        throw new InvalidOperationException("The specified storm network node feature class is invalid. It does not have an 'Id' field.");

      int progress = 0;
      int total = _nodesClass.FeatureCount(null);
      int increment = _CalcIncrement(total);

      IFeatureCursor cursor = null;
      try
      {
        cursor = _nodesClass.Search(null, false);
        IFeature node = cursor.NextFeature();
        while (node != null)
        {
          try
          {
            int id = (int)node.get_Value(idField);
            double load = _SafeConvertToDouble(node.get_Value(loadField), 0.0);
            if (_stormNetwork.Nodes.Contains(id))
            {
              _stormNetwork.Nodes[id].TotalAccumulatedLoad = load * conversionFactor;
            }
          }
          finally
          {
            ComUtility.ReleaseComObject(node);
          }

          if (++progress % increment == 0)
          {
            ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading accumulated load values...", true);
            OnProgress(eventArgs);
            if (eventArgs.Cancel)
              throw new UserCanceledException();
          }

          node = cursor.NextFeature();
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} load values read.", progress)), true);
    }

    public void LoadTotalAccumulation(string flowAccumulationField, string loadAccumulationField, double flowConversionFactor = 1.0, double loadConversionFactor = 1.0)
    {
      LoadTotalFlowAccumulation(flowAccumulationField, flowConversionFactor);
      LoadTotalLoadAccumulation(loadAccumulationField, loadConversionFactor);
    }

    public void UpdateStormNetwork(double flowConversionFactor = 1.0, double loadConversionFactor = 1.0)
    {
      if (_stormNetwork != null && _nodesClass != null && _conveyanceClass != null)
      {
        _UpdateStormNetworkNodes(_nodesClass, flowConversionFactor, loadConversionFactor);
        _UpdateStormNetworkConveyance(_conveyanceClass, flowConversionFactor, loadConversionFactor);
      }
    }

    private void _InitializeStormNetwork(StormDrainNetwork stormNetwork)
    {
      // Extract features
      Dictionary<int, Catchment> catchments = _ReadCatchments(true);
      _ReadInlets(stormNetwork, catchments.Keys.ToList());
      _ReadOutlets(stormNetwork);

      // Create links between network elements
      _LinkCatchmentsToInlets(catchments, stormNetwork.Inlets);
      _LinkNodesUnderground(stormNetwork);
      _LinkNodesOverland(stormNetwork);

      foreach (Catchment catchment in stormNetwork.Catchments)
        catchment.Polygon = null;
    }

    private Dictionary<int, Catchment> _ReadCatchments(bool includeShapes)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Reading catchments...", true), true);

      Dictionary<int, Catchment> catchments = new Dictionary<int, Catchment>();

      int progress = 0;
      int total = _catchmentClass.FeatureCount(null);
      int increment = _CalcIncrement(total);

      catchments.Clear();
      IFeatureCursor cursor = null;
      try
      {
        int gridCodeField = _catchmentClass.FindField("GridCode");
        cursor = _catchmentClass.Search(null, false);
        IFeature feature = cursor.NextFeature();
        while (feature != null)
        {
          int gridCode = Convert.ToInt32(feature.get_Value(gridCodeField));

          Catchment catchment;
          if (catchments.TryGetValue(gridCode, out catchment))
          {
            if (includeShapes)
              catchment.Polygon = ((ITopologicalOperator)catchment.Polygon).Union(feature.ShapeCopy) as IPolygon;
          }
          else
          {
            catchment = new Catchment();
            if (includeShapes)
              catchment.Polygon = feature.ShapeCopy as IPolygon;
            catchments.Add(gridCode, catchment);
          }

          if (++progress % increment == 0)
          {
            ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading catchments...", true);
            OnProgress(eventArgs);
            if (eventArgs.Cancel)
              throw new UserCanceledException();
          }

          feature = cursor.NextFeature();
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} catchments read.", catchments.Count)), true);

      return catchments;
    }

    private void _ReadInlets(StormDrainNetwork stormNetwork, List<int> inletIDs)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Reading inlets...", true), true);

      INetwork network = _geometricNetwork.Network;
      INetElements netElements = network as INetElements;

      int progress = 0;
      int total = inletIDs.Count;
      int increment = _CalcIncrement(total);

      int globalId = -1;
      int globalIdField = -1;
      object dbValue = null;

      int added = 0;
      foreach (int inletID in inletIDs)
      {
        int classID, userID, subID;
        netElements.QueryIDs(inletID, esriElementType.esriETJunction, out classID, out userID, out subID);

        INetElementDescriptionEdit nodeDescription = new NetElementDescriptionClass();
        nodeDescription.ElementType_2 = esriElementType.esriETJunction;
        nodeDescription.UserClassID_2 = classID;
        nodeDescription.UserID_2 = userID;
        nodeDescription.UserSubID_2 = subID;

        INetworkFeature netFeature = _geometricNetwork.get_NetworkFeature(nodeDescription);
        IFeature node = netFeature as IFeature;
        
        if (node != null)
        {
          try
          {
            string className = ((IDataset)node.Class).Name;

            // Reusing instance of class name reduces memory footprint slightly (yes, this works...)
            string existing;
            if (_classNames.TryGetValue(className, out existing))
            {
              // This will return false, meaning new instance of class name with same value as existing
              //System.Diagnostics.Debug.WriteLine(object.ReferenceEquals(className, existing));

              // Thus, this reduces memory when used in FeatureIdentifiers below since they all point
              // to the same string instance for the same feature class. Also the original className 
              // retrieved above from the feature class will be gargbage collected.
              className = existing;
            }
            else
            {
              _classNames.Add(className, className);
            }

            globalId = -1;
            globalIdField = node.Fields.FindField(GLOBALID_FIELD);
            if (globalIdField != -1)
            {
              dbValue = node.get_Value(globalIdField);
              globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
            }

            InletNode inletNode = new InletNode(inletID);
            inletNode.DataSource = new FeatureIdentifier() { FeatureClassName = className, ObjectId = userID, GlobalId = globalId };
            inletNode.Point = node.ShapeCopy as IPoint;
            inletNode.Enabled = !_IsDisabled(netFeature);

            stormNetwork.Nodes.Add(inletNode);
            added++;

            if (++progress % increment == 0)
            {
              ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading inlets...", true);
              OnProgress(eventArgs);
              if (eventArgs.Cancel)
                throw new UserCanceledException();
            }
          }
          finally
          {
            ComUtility.ReleaseComObject(node);
          }
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} inlets read.", added)), true);
    }

    private void _ReadOutlets(StormDrainNetwork stormNetwork)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Finding outlets...", true), true);

      INetwork network = _geometricNetwork.Network;
      INetTopology netTopology = network as INetTopology;

      int progress = 0;
      int total = 0;

      IEnumFeatureClass junctionClasses = _geometricNetwork.get_ClassesByType(esriFeatureType.esriFTSimpleJunction);
      IFeatureClass junctionClass = junctionClasses.Next();
      while (junctionClass != null)
      {
        total += junctionClass.FeatureCount(null);
        ComUtility.ReleaseComObject(junctionClass);
        junctionClass = junctionClasses.Next();
      }

      int increment = _CalcIncrement(total);

      int added = 0;
      junctionClasses.Reset();
      junctionClass = junctionClasses.Next();
      while (junctionClass != null)
      {
        IFeatureCursor cursor = null;
        try
        {
          int edgeId;
          bool towardJunction;
          
          string className = ((IDataset)junctionClass).Name;

          int globalId = -1;
          int globalIdField = junctionClass.FindField(GLOBALID_FIELD);
          object dbValue = null;

          cursor = junctionClass.Search(null, false);
          IFeature node = cursor.NextFeature();
          while (node != null)
          {
            try
            {
              ISimpleJunctionFeature junction = node as ISimpleJunctionFeature;
              bool isDownstream = junction.EdgeFeatureCount > 0;
              for (int i = 0; i < junction.EdgeFeatureCount; i++)
              {
                netTopology.GetAdjacentEdge(junction.EID, i, out edgeId, out towardJunction);
                if (!towardJunction)
                {
                  isDownstream = false;
                  break;
                }
              }

              if (isDownstream)
              {
                globalId = -1;
                if (globalIdField != -1)
                {
                  dbValue = node.get_Value(globalIdField);
                  globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
                }

                OutletNode outletNode = new OutletNode(junction.EID);
                outletNode.DataSource = new FeatureIdentifier() { FeatureClassName = className, ObjectId = node.OID, GlobalId = globalId};
                outletNode.Point = node.ShapeCopy as IPoint;
                outletNode.Enabled = !_IsDisabled(node as INetworkFeature);

                stormNetwork.Nodes.Add(outletNode);
                added++;
              }

              if (++progress % increment == 0)
              {
                ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Finding outlets...", true);
                OnProgress(eventArgs);
                if (eventArgs.Cancel)
                  throw new UserCanceledException();
              }
            }
            finally
            {
              ComUtility.ReleaseComObject(node);
            }
            node = cursor.NextFeature();
          }
        }
        finally
        {
          ComUtility.ReleaseComObject(cursor);
        }

        ComUtility.ReleaseComObject(junctionClass);
        junctionClass = junctionClasses.Next();
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} outlets found.", added)), true);
    }

    private void _LinkCatchmentsToInlets(Dictionary<int, Catchment> catchments, ElementCollection<InletNode> inletNodes)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Linking catchments to inlets..."), true);

      int progress = 0;
      int total = inletNodes.Count;
      int increment = _CalcIncrement(total);

      foreach (int inletID in inletNodes.IdList)
      {
        inletNodes[inletID].CatchmentArea = catchments[inletID];

        if (++progress % increment == 0)
        {
          OnProgress(new ProgressEventArgs(progress, total, "Linking catchments to inlets..."));
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, "Catchments and inlets linked."), true);
    }

    private void _LinkNodesUnderground(StormDrainNetwork stormNetwork)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Reading pipe shapes...", true), true);

      INetElements netElements = _geometricNetwork.Network as INetElements;
      IForwardStarGEN forwardStar = _geometricNetwork.Network.CreateForwardStar(false, null, null, null, null) as IForwardStarGEN;

      int progress = 0;
      int total = 0;

      IEnumFeatureClass edgeClasses = _geometricNetwork.get_ClassesByType(esriFeatureType.esriFTSimpleEdge);
      IFeatureClass edgeClass = edgeClasses.Next();
      while (edgeClass != null)
      {
        try
        {
          total += edgeClass.FeatureCount(null);
        }
        finally
        {
          ComUtility.ReleaseComObject(edgeClass);
        }
        edgeClass = edgeClasses.Next();
      }

      edgeClasses = _geometricNetwork.get_ClassesByType(esriFeatureType.esriFTComplexEdge);
      edgeClass = edgeClasses.Next();
      while (edgeClass != null)
      {
        try
        {
          total += edgeClass.FeatureCount(null);
        }
        finally
        {
          ComUtility.ReleaseComObject(edgeClass);
        }
        edgeClass = edgeClasses.Next();
      }

      int increment = _CalcIncrement(total);

      // Cache shapes for performance
      Dictionary<int, IPolyline> lineDictionary = new Dictionary<int, IPolyline>(total);
      Dictionary<int, FeatureIdentifier> idDictionary = new Dictionary<int, FeatureIdentifier>(total);

      edgeClasses = _geometricNetwork.get_ClassesByType(esriFeatureType.esriFTSimpleEdge);
      edgeClass = edgeClasses.Next();
      while (edgeClass != null)
      {
        try
        {
          string edgeClassName = ((IDataset)edgeClass).Name;

          int globalIdField = edgeClass.FindField(GLOBALID_FIELD);
          int globalId = -1;
          object dbValue = null;

          IFeatureCursor cursor = edgeClass.Search(null, true);
          try
          {
            IFeature edgeFeature = cursor.NextFeature();
            while (edgeFeature != null)
            {
              try
              {
                int eid = ((ISimpleEdgeFeature)edgeFeature).EID;
                lineDictionary.Add(eid, (IPolyline)edgeFeature.ShapeCopy);

                globalId = -1;
                if (globalIdField != -1)
                {
                  dbValue = edgeFeature.get_Value(globalIdField);
                  globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
                }

                idDictionary.Add(eid, new FeatureIdentifier() { FeatureClassName = edgeClassName, ObjectId = edgeFeature.OID , GlobalId = globalId});
              }
              finally
              {
                ComUtility.ReleaseComObject(edgeFeature);
              }

              if (++progress % increment == 0)
              {
                ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading pipe shapes...", true);
                OnProgress(eventArgs);
                if (eventArgs.Cancel)
                  throw new UserCanceledException();
              }

              edgeFeature = cursor.NextFeature();
            }
          }
          finally
          {
            ComUtility.ReleaseComObject(cursor);
          }
        }
        finally
        {
          ComUtility.ReleaseComObject(edgeClass);
        }
        edgeClass = edgeClasses.Next();
      }

      edgeClasses = _geometricNetwork.get_ClassesByType(esriFeatureType.esriFTComplexEdge);
      edgeClass = edgeClasses.Next();
      while (edgeClass != null)
      {
        try
        {
          string edgeClassName = ((IDataset)edgeClass).Name;

          int globalIdField = edgeClass.FindField(GLOBALID_FIELD);
          int globalId = -1;
          object dbValue = null;

          IFeatureCursor cursor = edgeClass.Search(null, true);
          try
          {
            IFeature edgeFeature = cursor.NextFeature();
            while (edgeFeature != null)
            {
              try
              {
                globalId = -1;
                if (globalIdField != -1)
                {
                  dbValue = edgeFeature.get_Value(globalIdField);
                  globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
                }

                FeatureIdentifier featureId = new FeatureIdentifier() { FeatureClassName = edgeClassName, ObjectId = edgeFeature.OID, GlobalId = globalId };
                IEnumNetEID eids = netElements.GetEIDs(edgeClass.ObjectClassID, edgeFeature.OID, esriElementType.esriETEdge);
                int eid = eids.Next();
                while (eid != 0)
                {
                  IPolyline lineGeometry = null;
                  try
                  {
                    lineGeometry = (IPolyline)((IComplexEdgeFeature)edgeFeature).get_GeometryForEID(eid);
                  }
                  catch (System.Runtime.InteropServices.COMException comex)
                  {
                    string message = string.Format("Unable to get geometry for EID {0} (FC: {1} OID: {2}). Using straight line instead.", eid, edgeClassName, edgeFeature.OID);
                    OnWarning(new WarningEventArgs(message, comex));
                    lineGeometry = new PolylineClass();
                  }

                  lineDictionary.Add(eid, lineGeometry);
                  idDictionary.Add(eid, featureId);
                  eid = eids.Next();
                }
              }
              finally
              {
                ComUtility.ReleaseComObject(edgeFeature);
              }

              if (++progress % increment == 0)
              {
                ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading pipe shapes...", true);
                OnProgress(eventArgs);
                if (eventArgs.Cancel)
                  throw new UserCanceledException();
              }

              edgeFeature = cursor.NextFeature();
            }
          }
          finally
          {
            ComUtility.ReleaseComObject(cursor);
          }
        }
        finally
        {
          ComUtility.ReleaseComObject(edgeClass);
        }
        edgeClass = edgeClasses.Next();
      }

      OnProgress(new ProgressEventArgs(0, 1, "Generating underground links...", true), true);

      progress = 0;
      total = stormNetwork.Outlets.Count;
      increment = _CalcIncrement(total);

      stormNetwork.ResetTraversedNodes();
      foreach (OutletNode outletNode in stormNetwork.Outlets)
      {
        _AddUpstreamElements(outletNode, forwardStar, netElements, lineDictionary, idDictionary);

        if (++progress % increment == 0)
        {
          ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Generating underground links...", true);
          OnProgress(eventArgs);
          if (eventArgs.Cancel)
            throw new UserCanceledException();
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, "Underground links generated."), true);
    }

    private void _AddUpstreamElements(NodePoint node, IForwardStarGEN forwardStar, INetElements netElements, Dictionary<int, IPolyline> lineDictionary, Dictionary<int, FeatureIdentifier> idDictionary)
    {
      // Break out of loops
      if (node.Traversed)
        return;
      else
        node.Traversed = true;

      int edgeCount = 0;
      forwardStar.FindAdjacent(0, node.ID, out edgeCount);

      if (edgeCount > 0)
      {
        int[] edgeEIDs = new int[edgeCount];
        int[] junctionEIDs = new int[edgeCount];
        bool[] towardsJunction = new bool[edgeCount];
        object[] weights = new object[edgeCount];

        forwardStar.QueryAdjacentEdges(ref edgeEIDs, ref towardsJunction, ref weights);
        forwardStar.QueryAdjacentJunctions(ref junctionEIDs, ref weights);
        for (int i = 0; i < edgeCount; i++)
        {
          if (towardsJunction[i])
          {
            bool edgeEnabled = !_IsDisabled(edgeEIDs[i], esriElementType.esriETEdge, netElements);

            if (edgeEnabled)
            {
              Conveyance conveyance = new Conveyance(edgeEIDs[i]);
              conveyance.DataSource = idDictionary[edgeEIDs[i]];
              conveyance.Line = lineDictionary[edgeEIDs[i]];
              conveyance.Enabled = edgeEnabled;
              conveyance.IsSubSurface = true;
              conveyance.DownstreamNode = node;

              NodePoint upNode = null;
              if (node.Network.Nodes.Contains(junctionEIDs[i]))
              {
                upNode = node.Network.Nodes[junctionEIDs[i]];
              }
              else
              {
                FeatureIdentifier source = null;
                IFeature junctionFeature = (IFeature)_GetNetworkFeature(junctionEIDs[i], esriElementType.esriETJunction, netElements);
                if (junctionFeature != null)
                {
                  string className = ((IDataset)junctionFeature.Class).Name;

                  // Reusing instance of class name reduces memory footprint slightly (yes, this works...)
                  string existing;
                  if (_classNames.TryGetValue(className, out existing))
                  {
                    // This will return false, meaning new instance of class name with same value as existing
                    //System.Diagnostics.Debug.WriteLine(object.ReferenceEquals(className, existing));

                    // Thus, this reduces memory when used in FeatureIdentifiers below since they all point
                    // to the same string instance for the same feature class. Also the original className 
                    // retrieved above from the feature class will be gargbage collected.
                    className = existing;
                  }
                  else
                  {
                    _classNames.Add(className, className);
                  }

                  int globalId = -1;
                  int globalIdField = junctionFeature.Fields.FindField(GLOBALID_FIELD);
                  if (globalIdField != -1)
                  {
                    object dbValue = junctionFeature.get_Value(globalIdField);
                    globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
                  }

                  source = new FeatureIdentifier() { FeatureClassName = className, ObjectId = junctionFeature.OID, GlobalId = globalId };
                }

                upNode = new NodePoint(junctionEIDs[i]);
                upNode.DataSource = source;
                upNode.Point = _geometricNetwork.get_GeometryForJunctionEID(junctionEIDs[i]) as IPoint;
                upNode.Enabled = !_IsDisabled(junctionEIDs[i], esriElementType.esriETJunction, netElements);
                node.Network.Nodes.Add(upNode);
              }
              conveyance.UpstreamNode = upNode;

              // For cases where geometry could not be read from geometric network due to an error.
              if (conveyance.Line == null || conveyance.Line.IsEmpty)
              {
                conveyance.Line = new PolylineClass();
                conveyance.Line.FromPoint = conveyance.UpstreamNode.Point;
                conveyance.Line.ToPoint = conveyance.DownstreamNode.Point;
              }

              node.Network.Conveyances.Add(conveyance);

              _AddUpstreamElements(conveyance.UpstreamNode, forwardStar, netElements, lineDictionary, idDictionary);
            }
          }
        }
      }
    }

    private void _LinkNodesOverland(StormDrainNetwork stormNetwork)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Generating overland links...", true), true);

      Dictionary<long, List<OutletNode>> indexedOutlets = new Dictionary<long, List<OutletNode>>();
      ElementCollection<OutletNode> remainingOutlets = new ElementCollection<OutletNode>();
      foreach (OutletNode outlet in stormNetwork.Outlets)
      {
        if (outlet.Enabled)
        {
          long index = -1;
          try
          {
            index = _GetRasterIndex(outlet.Point);
          }
          catch (IndexOutOfRangeException)
          {
            string message = string.Format("{0} is outside the flow direction grid. No downstream overland link created.", outlet);
            OnWarning(new WarningEventArgs(message));
          }

          // Only include outlets with valid flow direction
          if (index != -1 && !double.IsNaN(RasterUtility.GetPixelValueDouble(_flowDirRaster, outlet.Point)))
          {
            if (indexedOutlets.ContainsKey(index))
              indexedOutlets[index].Add(outlet);
            else
              indexedOutlets[index] = new List<OutletNode>(new OutletNode[] { outlet });

            remainingOutlets.Add(outlet);
          }
        }
      }

      int progress = 0;
      int total = remainingOutlets.Count;
      int increment = _CalcIncrement(total);

      long w = ((IRasterProps)_flowDirRaster).Width;

      foreach (InletNode inletNode in stormNetwork.Inlets)
      {
        int startRow, startColumn, width, height;
        _GetCatchmentCellIndeces(inletNode.CatchmentArea.Polygon.Envelope, out startRow, out startColumn, out width, out height);
        IRelationalOperator catchmentShape = inletNode.CatchmentArea.Polygon as IRelationalOperator;

        for (long r = startRow; r < startRow + height; r++)
        {
          for (long c = startColumn; c < startColumn + width; c++)
          {
            long i = r * w + c;

            if (indexedOutlets.ContainsKey(i))
            {
              foreach (OutletNode outletNode in indexedOutlets[i])
              {
                if (catchmentShape.Contains(outletNode.Point))
                {
                  IPolyline line = _GetFlowLineForOutlet(outletNode, inletNode.CatchmentArea.Polygon.Envelope);
                  if (line == null)
                  {
                    line = new PolylineClass();
                    line.FromPoint = outletNode.Point;
                    line.ToPoint = inletNode.Point;
                  }

                  // COMMENTED OUT: Snapping can lead to inaccurate conversion to raster later (causing zero-length conveyances)
                  // snap ends to pipe network
                  //line.FromPoint = outletNode.Point;
                  //line.ToPoint = inletNode.Point;

                  FlowLine flowLine = new FlowLine();
                  flowLine.DownstreamNode = inletNode;
                  flowLine.UpstreamNode = outletNode;
                  flowLine.Line = line;

                  stormNetwork.Conveyances.Add(flowLine);

                  remainingOutlets.Remove(outletNode);

                  if (++progress % increment == 0)
                  {
                    ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Generating overland links...", true);
                    OnProgress(eventArgs);
                    if (eventArgs.Cancel)
                      throw new UserCanceledException();
                  }
                }
              }
            }
          }
        }
      }

      // Process remaining outlets (these are terminal outlets--not associated with existing downstream inlets)
      OnProgress(new ProgressEventArgs(progress, total, "Generating overland links (terminal)...", false), true);
      if (_extractTerminalFlow)
      {
        OutletNode[] terminalOutlets = remainingOutlets.ToArray();
        IPolyline[] lines = _GetFlowLinesForOutlets(terminalOutlets);
        for (int i = 0; i < lines.Length; i++)
        {
          if (lines[i] != null && !lines[i].IsEmpty)
          {
            // COMMENTED OUT: Snapping can lead to inaccurate conversion to raster later (causing zero-length conveyances)
            // snap end to pipe network
            //lines[i].FromPoint = terminalOutlets[i].Point;

            FlowLine flowLine = new FlowLine();
            flowLine.UpstreamNode = terminalOutlets[i];
            flowLine.Line = lines[i];

            stormNetwork.Conveyances.Add(flowLine);
          }
          else
          {
            Logger.Log(string.Format("StormDrainNetworkLoader._LinkNodesOverland: {0} does not have flowline.", terminalOutlets[i]));
          }

          if (++progress % increment == 0)
          {
            ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Generating overland links (terminal)...", true);
            OnProgress(eventArgs);
            if (eventArgs.Cancel)
              throw new UserCanceledException();
          }
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, "Overland links generated."), true);
    }

    private void _GetCatchmentCellIndeces(IEnvelope envelope, out int startRow, out int startColumn, out int width, out int height)
    {
      double cellSize = ((IRasterAnalysisProps)_flowDirRaster).PixelWidth;

      IPoint origin = envelope.UpperLeft;
      origin.X += 0.5 * cellSize;
      origin.Y -= 0.5 * cellSize;

      ((IRaster2)_flowDirRaster).MapToPixel(origin.X, origin.Y, out startColumn, out startRow);

      IPoint right = envelope.UpperRight;
      right.X -= 0.5 * cellSize;
      right.Y -= 0.5 * cellSize;

      int colRight, rowRight;
      ((IRaster2)_flowDirRaster).MapToPixel(right.X, right.Y, out colRight, out rowRight);

      IPoint bottom = envelope.LowerLeft;
      bottom.X += 0.5 * cellSize;
      bottom.Y += 0.5 * cellSize;

      int colBottom, rowBottom;
      ((IRaster2)_flowDirRaster).MapToPixel(bottom.X, bottom.Y, out colBottom, out rowBottom);

      width = (colRight - startColumn) + 1;
      height = (rowBottom - startRow) + 1; //rows indexed from top
    }

    private IPolyline _GetFlowLineForOutlet(OutletNode outlet, IEnvelope extentHint = null)
    {
      // Get flow direction
      IGeoDataset flowDirDataset = (IGeoDataset)_flowDirRaster;

      // Check if outlet point is within flowdir extent. If not, just return null.
      if (!_IsWithin(outlet.Point, flowDirDataset.Extent))
        return null;

      // Get flow path geometry
      IPolyline flowPath = null;
      IDistanceOp distanceOp = new RasterDistanceOpClass();
      string outputDir = null;
      try
      {
        if (extentHint == null && outlet.ReceivingCatchment != null)
        {
          extentHint = outlet.ReceivingCatchment.Polygon.Envelope;
        }

        if (extentHint == null)
        {
          if (_extractTerminalFlow)
            this.SetAnalysisEnvironment((IRasterAnalysisEnvironment)distanceOp, _flowDirRaster);
          else
            return null;
        }
        else
        {
          // Expand to ensure full boundary of catchment is included.
          double cellWidth = ((IRasterAnalysisProps)_flowDirRaster).PixelWidth;
          double cellHeight = ((IRasterAnalysisProps)_flowDirRaster).PixelHeight;
          extentHint.Expand(cellWidth * 3.0, cellHeight * 3.0, false);

          this.SetAnalysisEnvironment((IRasterAnalysisEnvironment)distanceOp, extentHint);
        }
        outputDir = ((IRasterAnalysisEnvironment)distanceOp).OutWorkspace.PathName;

        // Process point
        IPointCollection pointCollection = new MultipointClass();
        ((IGeometry)pointCollection).SpatialReference = flowDirDataset.SpatialReference;
        IPoint outletPoint = ((IClone)outlet.Point).Clone() as IPoint;
        outletPoint.Project(flowDirDataset.SpatialReference);
        pointCollection.AddPoint(outletPoint);

        IGeometryCollection flowPaths = distanceOp.CostPathAsPolyline(pointCollection, flowDirDataset, flowDirDataset);
        if (flowPaths.GeometryCount > 0)
          flowPath = flowPaths.get_Geometry(0) as IPolyline;
      }
      catch
      {
        Logger.Log(string.Format("Failed to get flow path for outlet {0}.", outlet.ID));
        throw;
      }
      finally
      {
        ComUtility.ReleaseComObject(distanceOp);
        try
        {
          if (Directory.Exists(outputDir))
            Directory.Delete(outputDir, true);
        }
        catch (Exception ex)
        {
          Logger.Log("StormDrainNetworkLoader._GetFlowLineForOutlet: (Warning) Cannot delete temp directory: " + ex.GetType().FullName + ": " + ex.Message);
        }
      }

      return flowPath;
    }

    private IPolyline[] _GetFlowLinesForOutlets(OutletNode[] outlets)
    {
      // Get flow direction
      IGeoDataset flowDirDataset = (IGeoDataset)_flowDirRaster;

      // Check if outlet point is within flowdir extent. If not, just return null.
      //if (!_IsWithin(outlet.Point, flowDirDataset.Extent))
      //  return null;

      // Get flow path geometry
      IPolyline[] flowPaths = new IPolyline[outlets.Length];
      IDistanceOp distanceOp = new RasterDistanceOpClass();
      string outputDir = null;
      try
      {
        this.SetAnalysisEnvironment((IRasterAnalysisEnvironment)distanceOp, _flowDirRaster);
        outputDir = ((IRasterAnalysisEnvironment)distanceOp).OutWorkspace.PathName;

        // Process points
        IPointCollection pointCollection = new MultipointClass();
        ((IGeometry)pointCollection).SpatialReference = flowDirDataset.SpatialReference;
        for (int i = 0; i < outlets.Length; i++)
        {
          IPoint outletPoint = ((IClone)outlets[i].Point).Clone() as IPoint;
          outletPoint.Project(flowDirDataset.SpatialReference);
          pointCollection.AddPoint(outletPoint);
        }

        IGeometryCollection flowPathCollection = distanceOp.CostPathAsPolyline(pointCollection, flowDirDataset, flowDirDataset);
        if (flowPathCollection.GeometryCount != outlets.Length)
          Logger.Log("StormDrainNetworkLoader._GetFlowLinesForOutlets: Resulting terminal flowline count does not match terminal outlet count.");

        for (int i = 0; i < outlets.Length; i++)
        {
          flowPaths[i] = flowPathCollection.get_Geometry(i) as IPolyline;
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(distanceOp);
        try
        {
          if (Directory.Exists(outputDir))
            Directory.Delete(outputDir, true);
        }
        catch (Exception ex)
        {
          Logger.Log("StormDrainNetworkLoader._GetFlowLinesForOutlets: (Warning) Cannot delete temp directory: " + ex.GetType().FullName + ": " + ex.Message);
        }
      }

      return flowPaths;
    }

    private static bool _IsWithin(IPoint point, IEnvelope envelope)
    {
      return (point.X >= envelope.LowerLeft.X &&
          point.X <= envelope.LowerRight.X &&
          point.Y >= envelope.LowerLeft.Y &&
          point.Y <= envelope.UpperLeft.Y);
    }

    private IRaster _CalculateLocalFlowAccumulation()
    {
      OnProgress(new ProgressEventArgs(0, 1, "Calculating local flow accumulation..."), true);

      Geoprocessor gp = GeoprocessingUtility.GetGeoprocessor(true, false, true, _flowDirRaster);
      try
      {
        IDataset flowDirDataset = (IDataset)((IRasterAnalysisProps)_flowDirRaster).RasterDataset;
        string flowDirPath = System.IO.Path.Combine(flowDirDataset.Workspace.PathName, flowDirDataset.Name);

        IDataset runoffDataset = (IDataset)((IRasterAnalysisProps)_runoffRaster).RasterDataset;
        string runoffPath = System.IO.Path.Combine(runoffDataset.Workspace.PathName, runoffDataset.Name);

        FlowAccumulation flowAccTool = new FlowAccumulation();
        flowAccTool.in_flow_direction_raster = flowDirPath;
        flowAccTool.in_weight_raster = runoffPath;
        flowAccTool.out_accumulation_raster = "local_runoff_accumulation";

        object result = GeoprocessingUtility.RunGpTool(gp, flowAccTool, OnProgress);
        Logger.Log(string.Format("Local flow accumulation: {0}", result));

        OnProgress(new ProgressEventArgs(0, 1, "Returning the local runoff accumulation raster..."), true);

        string workspacePath = System.IO.Path.GetDirectoryName((string)result);
        string rasterName = System.IO.Path.GetFileName((string)result);

        IWorkspaceFactory wsf = new FileGDBWorkspaceFactoryClass();
        IWorkspace ws = wsf.OpenFromFile(workspacePath, 0);
        IRasterDataset rasterDataset = ((IRasterWorkspaceEx)ws).OpenRasterDataset(rasterName);

        IRaster raster = rasterDataset.CreateDefaultRaster();
        return raster;
      }
      finally
      {
        GeoprocessingUtility.ResetGeoprocessor();
      }
    }

    private void _ExtractLocalAccumulation(ICollection<Catchment> catchments, IRaster localAccumulationRaster)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Extracting local flow accumulation...", true), true);

      int progress = 0;
      int total = catchments.Count;
      int increment = _CalcIncrement(total);

      // order for performance (not a huge gain, unfortunately...)
      List<Catchment> orderedCatchments = new List<Catchment>(catchments);
      orderedCatchments.Sort((c1, c2) => c2.DrainageNode.Point.Y.CompareTo(c1.DrainageNode.Point.Y));

      foreach (Catchment catchment in orderedCatchments)
      {
        double localFlowAcc = RasterUtility.GetPixelValueDouble(localAccumulationRaster, catchment.DrainageNode.Point);
        catchment.LocalFlowAccumulation = double.IsNaN(localFlowAcc) ? 0.0 : localFlowAcc;

        if (++progress % increment == 0)
        {
          ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Extracting local flow accumulation...", true);
          OnProgress(eventArgs);
          if (eventArgs.Cancel)
            throw new UserCanceledException();
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, "Local flow accumulation extracted."), true);
    }

    private void _SaveStormNetwork(IWorkspace workspace, string baseName, double flowConversionFactor = 1.0)
    {
      IFeatureClass nodeClass = _MakeNodeFeatureClass(workspace, baseName + "_nodes");
      _UpdateStormNetworkNodes(nodeClass, flowConversionFactor);
      _nodesClass = nodeClass;

      IFeatureClass conveyanceClass = _MakeConveyanceFeatureClass(workspace, baseName + "_conveyance");
      _UpdateStormNetworkConveyance(conveyanceClass, flowConversionFactor);
      _conveyanceClass = conveyanceClass;
    }

    private IFeatureClass _MakeNodeFeatureClass(IWorkspace workspace, string name)
    {
      ISpatialReference spatialRef = _GetPreferedSpatialReference();

      IFeatureClass nodeClass = FeatureClassUtility.CreateFeatureClass(workspace, name, esriGeometryType.esriGeometryPoint, spatialRef, ID_FIELD, true);
      FeatureClassUtility.AddField(ref nodeClass, TYPE_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref nodeClass, DOWNNODE_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref nodeClass, UPNODES_FIELD, esriFieldType.esriFieldTypeString, 512, false);
      FeatureClassUtility.AddField(ref nodeClass, ENABLED_FIELD, esriFieldType.esriFieldTypeSmallInteger, false);
      FeatureClassUtility.AddField(ref nodeClass, ERRORS_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref nodeClass, SOURCECLASS_FIELD, esriFieldType.esriFieldTypeString, 80, false);
      FeatureClassUtility.AddField(ref nodeClass, SOURCEOID_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref nodeClass, GLOBALID_FIELD, esriFieldType.esriFieldTypeInteger, false);

      return nodeClass;
    }

    private IFeatureClass _MakeConveyanceFeatureClass(IWorkspace workspace, string name)
    {
      ISpatialReference spatialRef = _GetPreferedSpatialReference();

      IFeatureClass conveyanceClass = FeatureClassUtility.CreateFeatureClass(workspace, name, esriGeometryType.esriGeometryPolyline, spatialRef, ID_FIELD, true);
      FeatureClassUtility.AddField(ref conveyanceClass, UPNODE_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref conveyanceClass, DOWNNODE_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref conveyanceClass, ENABLED_FIELD, esriFieldType.esriFieldTypeSmallInteger, false);
      FeatureClassUtility.AddField(ref conveyanceClass, SUBSURFACE_FIELD, esriFieldType.esriFieldTypeSmallInteger, false);
      FeatureClassUtility.AddField(ref conveyanceClass, ERRORS_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref conveyanceClass, SOURCECLASS_FIELD, esriFieldType.esriFieldTypeString, 80, false);
      FeatureClassUtility.AddField(ref conveyanceClass, SOURCEOID_FIELD, esriFieldType.esriFieldTypeInteger, false);
      FeatureClassUtility.AddField(ref conveyanceClass, GLOBALID_FIELD, esriFieldType.esriFieldTypeInteger, false);
      return conveyanceClass;
    }

    private void _UpdateStormNetworkNodes(IFeatureClass nodeClass, double flowConversionFactor = 1.0, double loadConversionFactor = 1.0)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Saving storm network nodes..."), true);

      int progress = 0;
      int total = _stormNetwork.Nodes.Count;
      int increment = _CalcIncrement(total);

      int idField = nodeClass.FindField(ID_FIELD);
      int typeField = nodeClass.FindField(TYPE_FIELD);
      int downNodeField = nodeClass.FindField(DOWNNODE_FIELD);
      int upNodesField = nodeClass.FindField(UPNODES_FIELD);
      int enabledField = nodeClass.FindField(ENABLED_FIELD);
      int errorField = nodeClass.FindField(ERRORS_FIELD);
      int sourceClassField = nodeClass.FindField(SOURCECLASS_FIELD);
      int sourceOidField = nodeClass.FindField(SOURCEOID_FIELD);
      int globalIdField = nodeClass.FindField(GLOBALID_FIELD);
      int accumField = string.IsNullOrWhiteSpace(_flowAccumField) ? -1 : nodeClass.FindField(_flowAccumField);
      int loadField = string.IsNullOrWhiteSpace(_loadAccumField) ? -1 : nodeClass.FindField(_loadAccumField);
      int tagField = string.IsNullOrWhiteSpace(_tagField) ? -1 : nodeClass.FindField(_tagField);

      if (!string.IsNullOrWhiteSpace(_flowAccumField) && accumField == -1)
        accumField = FeatureClassUtility.AddField(ref nodeClass, _flowAccumField, esriFieldType.esriFieldTypeDouble, true);

      if (!string.IsNullOrWhiteSpace(_loadAccumField) && loadField == -1)
        loadField = FeatureClassUtility.AddField(ref nodeClass, _loadAccumField, esriFieldType.esriFieldTypeDouble, true);

      int upNodesFieldLength = nodeClass.Fields.get_Field(upNodesField).Length;

      HashSet<int> nodesToProcess = new HashSet<int>(_stormNetwork.Nodes.IdList);

      // First update existing nodes
      IFeatureCursor cursor = null;
      try
      {
        cursor = nodeClass.Update(null, false);
        IFeature nodeFeature = cursor.NextFeature();
        while (nodeFeature != null)
        {
          try
          {
            int id = (int)nodeFeature.get_Value(idField);
            if (nodesToProcess.Contains(id))
            {
              NodePoint node = _stormNetwork.Nodes[id];
              nodeFeature.Shape = node.Point;

              nodeFeature.set_Value(typeField, (int)node.Type);
              nodeFeature.set_Value(downNodeField, node.DownstreamNode == null ? DBNull.Value : (object)node.DownstreamNode.ID);

              string upNodes = _GetUpNodesAsString(node, upNodesFieldLength);
              nodeFeature.set_Value(upNodesField, upNodes == string.Empty ? DBNull.Value : (object)upNodes);

              nodeFeature.set_Value(enabledField, node.Enabled ? (short)1 : (short)0);
              nodeFeature.set_Value(errorField, (int)node.Errors);

              if (accumField != -1)
              {
                nodeFeature.set_Value(accumField, node.TotalAccumulatedFlow * flowConversionFactor);
              }

              if (loadField != -1)
              {
                nodeFeature.set_Value(loadField, node.TotalAccumulatedLoad * loadConversionFactor);
              }

              if (sourceClassField != -1 && sourceOidField != -1)
              {
                FeatureIdentifier source = node.DataSource as FeatureIdentifier;
                if (source != null)
                {
                  nodeFeature.set_Value(sourceClassField, source.FeatureClassName);
                  nodeFeature.set_Value(sourceOidField, source.ObjectId);
                }
              }

              if (globalIdField != -1)
              {
                FeatureIdentifier source = node.DataSource as FeatureIdentifier;
                if (source != null)
                {
                  nodeFeature.set_Value(globalIdField, source.GlobalId);
                }
              }

              if (tagField != -1)
              {
                nodeFeature.set_Value(tagField, node.Tag == null ? DBNull.Value : node.Tag);
              }

              cursor.UpdateFeature(nodeFeature);
              nodesToProcess.Remove(id);

              if (++progress % increment == 0)
              {
                OnProgress(new ProgressEventArgs(progress, total, "Saving storm network nodes..."));
              }
            }
            else
            {
              nodeFeature.Delete();
            }
          }
          catch (Exception ex)
          {
            Logger.Log("StormDrainNetworkLoader._UpdateStormNetworkNodes: Unable to update storm network node: " + ex.Message);
            throw;
          }
          finally
          {
            ComUtility.ReleaseComObject(nodeFeature);
          }

          nodeFeature = cursor.NextFeature();
        }
        cursor.Flush();
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      // Then insert new nodes
      if (nodesToProcess.Count > 0)
      {
        try
        {
          cursor = nodeClass.Insert(false);
          IFeatureBuffer buffer = nodeClass.CreateFeatureBuffer();

          foreach (int id in nodesToProcess)
          {
            NodePoint node = _stormNetwork.Nodes[id];
            buffer.Shape = node.Point;

            buffer.set_Value(idField, node.ID);
            buffer.set_Value(typeField, (int)node.Type);
            buffer.set_Value(downNodeField, node.DownstreamNode == null ? DBNull.Value : (object)node.DownstreamNode.ID);

            string upNodes = _GetUpNodesAsString(node, upNodesFieldLength);
            buffer.set_Value(upNodesField, upNodes == string.Empty ? DBNull.Value : (object)upNodes);

            buffer.set_Value(enabledField, node.Enabled ? (short)1 : (short)0);
            buffer.set_Value(errorField, (int)node.Errors);

            if (accumField != -1)
            {
              buffer.set_Value(accumField, node.TotalAccumulatedFlow * flowConversionFactor);
            }

            if (loadField != -1)
            {
              buffer.set_Value(loadField, node.TotalAccumulatedLoad * loadConversionFactor);
            }

            if (sourceClassField != -1 && sourceOidField != -1)
            {
              FeatureIdentifier source = node.DataSource as FeatureIdentifier;
              if (source != null)
              {
                buffer.set_Value(sourceClassField, source.FeatureClassName);
                buffer.set_Value(sourceOidField, source.ObjectId);
              }
              else
              {
                buffer.set_Value(sourceClassField, DBNull.Value);
                buffer.set_Value(sourceOidField, DBNull.Value);
              }
            }

            if (globalIdField != -1)
            {
              FeatureIdentifier source = node.DataSource as FeatureIdentifier;
              if (source != null)
              {
                buffer.set_Value(globalIdField, source.GlobalId);
              }
              else
              {
                buffer.set_Value(globalIdField, DBNull.Value);
              }
            }

            if (tagField != -1)
            {
              buffer.set_Value(tagField, node.Tag == null ? DBNull.Value : node.Tag);
            }

            cursor.InsertFeature(buffer);

            if (++progress % increment == 0)
            {
              OnProgress(new ProgressEventArgs(progress, total, "Saving storm network nodes..."));
            }
          }

          cursor.Flush();
        }
        finally
        {
          ComUtility.ReleaseComObject(cursor);
        }
      }

      Logger.Log(string.Format("Saved/Updated Storm Network Nodes: {0}", GeoprocessingUtility.GetFeatureClassPath(nodeClass)));
      OnProgress(new ProgressEventArgs(1, 1, "Storm network nodes saved."), true);
    }

    private string _GetUpNodesAsString(NodePoint node, int upNodesFieldLength)
    {
      string upNodes = string.Empty;
      string separator = string.Empty;
      foreach (NodePoint upNode in node.UpstreamNodes)
      {
        string currentNodes = upNodes;

        upNodes += separator + upNode.ID.ToString();
        separator = ";";

        if (upNodes.Length > upNodesFieldLength)
        {
          upNodes = currentNodes;
          Logger.Log(string.Format("StormDrainNetworkLoader._GetUpNodesAsString: Upstream node ids too long for {0}.", node));
          break;
        }
      }
      return upNodes;
    }

    private void _UpdateStormNetworkConveyance(IFeatureClass conveyanceClass, double flowConversionFactor = 1.0, double loadConversionFactor = 1.0)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Saving storm network conveyance..."), true);

      int progress = 0;
      int total = _stormNetwork.Conveyances.Count;
      int increment = _CalcIncrement(total);

      int idField = conveyanceClass.FindField(ID_FIELD);
      int downNodeField = conveyanceClass.FindField(DOWNNODE_FIELD);
      int upNodeField = conveyanceClass.FindField(UPNODE_FIELD);
      int enabledField = conveyanceClass.FindField(ENABLED_FIELD);
      int subSurfaceField = conveyanceClass.FindField(SUBSURFACE_FIELD);
      int errorField = conveyanceClass.FindField(ERRORS_FIELD);
      int sourceClassField = conveyanceClass.FindField(SOURCECLASS_FIELD);
      int sourceOidField = conveyanceClass.FindField(SOURCEOID_FIELD);
      int globalIdField = conveyanceClass.FindField(GLOBALID_FIELD);
      int accumField = string.IsNullOrWhiteSpace(_flowAccumField) ? -1 : conveyanceClass.FindField(_flowAccumField);
      int loadField = string.IsNullOrWhiteSpace(_loadAccumField) ? -1 : conveyanceClass.FindField(_loadAccumField);
      int tagField = string.IsNullOrEmpty(_tagField) ? -1 : conveyanceClass.FindField(_tagField);

      if (!string.IsNullOrWhiteSpace(_flowAccumField) && accumField == -1)
        accumField = FeatureClassUtility.AddField(ref conveyanceClass, _flowAccumField, esriFieldType.esriFieldTypeDouble, true);

      if (!string.IsNullOrWhiteSpace(_loadAccumField) && loadField == -1)
        loadField = FeatureClassUtility.AddField(ref conveyanceClass, _loadAccumField, esriFieldType.esriFieldTypeDouble, true);

      HashSet<int> conveyancesToProcess = new HashSet<int>(_stormNetwork.Conveyances.IdList);

      // First update existing conveyances
      IFeatureCursor cursor = null;
      try
      {
        cursor = conveyanceClass.Update(null, false);
        IFeature conveyanceFeature = cursor.NextFeature();
        while (conveyanceFeature != null)
        {
          try
          {
            int id = (int)conveyanceFeature.get_Value(idField);
            if (conveyancesToProcess.Contains(id))
            {
              Conveyance conveyance = _stormNetwork.Conveyances[id];
              conveyanceFeature.Shape = conveyance.Line;

              conveyanceFeature.set_Value(downNodeField, conveyance.DownstreamNode == null ? DBNull.Value : (object)conveyance.DownstreamNode.ID);
              conveyanceFeature.set_Value(upNodeField, conveyance.UpstreamNode == null ? DBNull.Value : (object)conveyance.UpstreamNode.ID);
              conveyanceFeature.set_Value(enabledField, conveyance.Enabled ? (short)1 : (short)0);
              conveyanceFeature.set_Value(subSurfaceField, conveyance.IsSubSurface ? (short)1 : (short)0);
              conveyanceFeature.set_Value(errorField, (int)conveyance.Errors);

              if (accumField != -1)
              {
                conveyanceFeature.set_Value(accumField, conveyance.MinAccumulatedFlow * flowConversionFactor);
              }

              if (loadField != -1)
              {
                conveyanceFeature.set_Value(loadField, conveyance.MinAccumulatedLoad * loadConversionFactor);
              }

              if (sourceClassField != -1 && sourceOidField != -1)
              {
                FeatureIdentifier source = conveyance.DataSource as FeatureIdentifier;
                if (source != null)
                {
                  conveyanceFeature.set_Value(sourceClassField, source.FeatureClassName);
                  conveyanceFeature.set_Value(sourceOidField, source.ObjectId);
                }
              }

              if (globalIdField != -1)
              {
                FeatureIdentifier source = conveyance.DataSource as FeatureIdentifier;
                if (source != null)
                {
                  conveyanceFeature.set_Value(globalIdField, source.GlobalId);
                }
              }

              if (tagField != -1)
              {
                conveyanceFeature.set_Value(tagField, conveyance.Tag == null ? DBNull.Value : conveyance.Tag);
              }

              cursor.UpdateFeature(conveyanceFeature);
              conveyancesToProcess.Remove(id);

              if (++progress % increment == 0)
              {
                OnProgress(new ProgressEventArgs(progress, total, "Saving storm network conveyance..."));
              }
            }
            else
            {
              conveyanceFeature.Delete();
            }
          }
          catch (Exception ex)
          {
            Logger.Log("StormDrainNetworkLoader._UpdateStormNetwork: Unable to update storm network conveyance: " + ex.Message);
            throw;
          }
          finally
          {
            ComUtility.ReleaseComObject(conveyanceFeature);
          }
          conveyanceFeature = cursor.NextFeature();
        }

        cursor.Flush();
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      // Then insert new conveyances
      if (conveyancesToProcess.Count > 0)
      {
        try
        {
          cursor = conveyanceClass.Insert(false);
          IFeatureBuffer buffer = conveyanceClass.CreateFeatureBuffer();

          foreach (int id in conveyancesToProcess)
          {
            Conveyance conveyance = _stormNetwork.Conveyances[id];
            buffer.Shape = conveyance.Line;

            buffer.set_Value(idField, conveyance.ID);
            buffer.set_Value(downNodeField, conveyance.DownstreamNode == null ? DBNull.Value : (object)conveyance.DownstreamNode.ID);
            buffer.set_Value(upNodeField, conveyance.UpstreamNode == null ? DBNull.Value : (object)conveyance.UpstreamNode.ID);
            buffer.set_Value(enabledField, conveyance.Enabled ? (short)1 : (short)0);
            buffer.set_Value(subSurfaceField, conveyance.IsSubSurface ? (short)1 : (short)0);
            buffer.set_Value(errorField, (int)conveyance.Errors);

            if (accumField != -1)
            {
              buffer.set_Value(accumField, conveyance.MinAccumulatedFlow * flowConversionFactor);
            }

            if (loadField != -1)
            {
              buffer.set_Value(loadField, conveyance.MinAccumulatedLoad * loadConversionFactor);
            }

            if (sourceClassField != -1 && sourceOidField != -1)
            {
              FeatureIdentifier source = conveyance.DataSource as FeatureIdentifier;
              if (source != null)
              {
                buffer.set_Value(sourceClassField, source.FeatureClassName);
                buffer.set_Value(sourceOidField, source.ObjectId);
              }
              else
              {
                buffer.set_Value(sourceClassField, DBNull.Value);
                buffer.set_Value(sourceOidField, DBNull.Value);
              }
            }

            if (globalIdField != -1)
            {
              FeatureIdentifier source = conveyance.DataSource as FeatureIdentifier;
              if (source != null)
              {
                buffer.set_Value(globalIdField, source.GlobalId);
              }
              else
              {
                buffer.set_Value(globalIdField, DBNull.Value);
              }
            }

            if (tagField != -1)
            {
              buffer.set_Value(tagField, conveyance.Tag == null ? DBNull.Value : conveyance.Tag);
            }

            cursor.InsertFeature(buffer);

            if (++progress % increment == 0)
            {
              OnProgress(new ProgressEventArgs(progress, total, "Saving storm network conveyance..."));
            }
          }

          cursor.Flush();
        }
        finally
        {
          ComUtility.ReleaseComObject(cursor);
        }
      }

      Logger.Log(string.Format("Saved/Updated Storm Network Conveyances: {0}", GeoprocessingUtility.GetFeatureClassPath(conveyanceClass)));
      OnProgress(new ProgressEventArgs(1, 1, "Storm network conveyance saved."), true);
    }

    public IFeatureClass _CalculateTotalFlowAccumulation(IRaster localAccumulation, int threshold)
    {
      // Write out Python script
      OnProgress(new ProgressEventArgs(0, 1, "Writing out Python script to calculate Total Flow Accumulation..."), true);

      string documents = Environment.GetFolderPath(Environment.SpecialFolder.MyDocuments);
      string tempDir = System.IO.Path.Combine(documents, "temp");
      string scriptPath = System.IO.Path.Combine(tempDir, "CalculateTotalAccumulation.py");
      Directory.CreateDirectory(tempDir);

      Assembly assembly = Assembly.GetExecutingAssembly();
      using (StreamReader templateReader = new StreamReader(assembly.GetManifestResourceStream("Bes.Swsp.ToolBox.Templates.CalculateTotalAccumulation.py")))
      {
        using (StreamWriter scriptWriter = new StreamWriter(scriptPath))
        {
          while (templateReader.Peek() != -1)
          {
            scriptWriter.WriteLine(templateReader.ReadLine());
          }
        }
      }

      // Execute script
      OnProgress(new ProgressEventArgs(0, 1, "Executing Python script to calculate Total Flow Accumulation..."), true);

      IDataset conveyanceDataset = (IDataset)this.StormNetworkConveyanceClass;
      string conveyancePath = System.IO.Path.Combine(conveyanceDataset.Workspace.PathName, conveyanceDataset.Name);

      IDataset localAccumDataset = (IDataset) ((IRasterAnalysisProps)localAccumulation).RasterDataset;
      string localAccumPath = System.IO.Path.Combine(localAccumDataset.Workspace.PathName, localAccumDataset.Name);

      IDataset flowDirDataset = (IDataset)((IRasterAnalysisProps)_flowDirRaster).RasterDataset;
      string flowDirPath = System.IO.Path.Combine(flowDirDataset.Workspace.PathName, flowDirDataset.Name);

      string scratchWorkspace = GeoprocessingUtility.ScratchWorkspace;
      if (string.IsNullOrEmpty(scratchWorkspace))
        scratchWorkspace = conveyanceDataset.Workspace.PathName;

      string logPath = Logger.LogFile;
      if (logPath == null)
        logPath = string.Empty;

      string arguments = string.Format("{0} {1} {2} {3} {4} {5}", conveyancePath.Quoted(), localAccumPath.Quoted(), flowDirPath.Quoted(), threshold.ToString(), scratchWorkspace.Quoted(), logPath.Quoted());

      Process pythonProcess = null;
      bool restartLog = false;
      if (Logger.IsLogging)
      {
        // Release log file to allow logging by script
        Logger.StopLogging();
        restartLog = true;
      }
      try
      {
        pythonProcess = PythonUtility.StartPythonScript(scriptPath, arguments);

        while (!pythonProcess.HasExited)
        {
          System.Threading.Thread.Sleep(500);
          OnProgress(new ProgressEventArgs(0, 1, "Executing Python script to calculate Total Flow Accumulation..."));
        }

        if (File.Exists(scriptPath))
          File.Delete(scriptPath);
      }
      finally
      {
        if (restartLog)
          Logger.StartLogging();
      }

      if (pythonProcess != null && pythonProcess.HasExited && pythonProcess.ExitCode != 0)
      {
        throw new Exception("Python returned an error while executing Total Flow Accumulation");
      }

      OnProgress(new ProgressEventArgs(0, 1, "Returning the stream line feature class..."), true);

      IWorkspaceFactory wsf = new FileGDBWorkspaceFactoryClass();
      IWorkspace ws = wsf.OpenFromFile(conveyanceDataset.Workspace.PathName, 0);
      IFeatureClass streamClass = ((IFeatureWorkspace)ws).OpenFeatureClass("accumulated_surface_lines");

      IDataset streamDataset = (IDataset)streamClass;
      string streamClassPath = System.IO.Path.Combine(streamDataset.Workspace.PathName, streamDataset.Name);
      OnProgress(new ProgressEventArgs(1, 1, "Stream feature class: " + streamClassPath), true);

      return streamClass;
    }

    public void _RebuildSurfaceNetwork(IFeatureClass streamClass)
    {
      StormDrainNetwork overlandNetwork = _ReadStreamLinks(streamClass);
      Dictionary<long, NodePoint> indexedSurfaceNodes = _IndexSurfaceNodes();
      Dictionary<NodePoint, NodePoint> nodeMapping = _SnapEndpoints(overlandNetwork, indexedSurfaceNodes);
      _ReplaceOverlandNetwork(overlandNetwork, nodeMapping);
    }

    private StormDrainNetwork _ReadStreamLinks(IFeatureClass streamClass)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Reading stream lines...", true), true);

      int progress = 0;
      int total = streamClass.FeatureCount(null);
      int increment = _CalcIncrement(total);
      bool cancel = false;

      StormDrainNetwork overlandNetwork = new StormDrainNetwork();

      int idField = streamClass.FindField("arcid");
      int fromField = streamClass.FindField("from_node");
      int toField = streamClass.FindField("to_node");
      int valueField = streamClass.FindField("grid_code");

      string className = ((IDataset)streamClass).Name;

      IFeatureCursor cursor = streamClass.Search(null, true);
      try
      {
        IFeature streamFeature = cursor.NextFeature();
        while (streamFeature != null)
        {
          try
          {
            int id = (int)streamFeature.get_Value(idField);
            int upNodeId = (int)streamFeature.get_Value(fromField);
            int downNodeId = (int)streamFeature.get_Value(toField);
            double value = _SafeConvertToDouble(streamFeature.get_Value(valueField));

            FlowLine flowLine = new FlowLine(id);
            flowLine.Line = streamFeature.ShapeCopy as IPolyline;
            flowLine.DataSource = new FeatureIdentifier() { FeatureClassName = className, ObjectId = streamFeature.OID, GlobalId = -1 };

            if (overlandNetwork.Nodes.Contains(upNodeId))
            {
              flowLine.UpstreamNode = overlandNetwork.Nodes[upNodeId];
            }
            else
            {
              flowLine.UpstreamNode = new NodePoint(upNodeId);
              flowLine.UpstreamNode.Point = flowLine.Line.FromPoint;
            }

            if (overlandNetwork.Nodes.Contains(downNodeId))
            {
              flowLine.DownstreamNode = overlandNetwork.Nodes[downNodeId];
            }
            else
            {
              flowLine.DownstreamNode = new NodePoint(downNodeId);
              flowLine.DownstreamNode.Point = flowLine.Line.ToPoint;
            }

            flowLine.UpstreamNode.TotalAccumulatedFlow = value;
            overlandNetwork.Conveyances.Add(flowLine);

            if (++progress % increment == 0)
            {
              ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Reading stream lines...", true);
              OnProgress(eventArgs);
              if (eventArgs.Cancel)
              {
                cancel = true;
                break;
              }
            }
          }
          finally
          {
            ComUtility.ReleaseComObject(streamFeature);
          }

          streamFeature = cursor.NextFeature();
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      OnProgress(cancel
        ? new ProgressEventArgs(0, 0, "Canceled. Changes rolled back.")
        : new ProgressEventArgs(1, 1, "Stream lines read."), true);

      return overlandNetwork;
    }

    private Dictionary<long, NodePoint> _IndexSurfaceNodes()
    {
      OnProgress(new ProgressEventArgs(0, 1, "Indexing nodes..."), true);

      Dictionary<long, NodePoint> indexedSurfaceNodes = new Dictionary<long, NodePoint>(_stormNetwork.Nodes.Count);
      foreach (NodePoint node in _stormNetwork.Nodes)
      {
        // Node is at surface if it is an inlet node (always), or if there is an overland conveyance linked to it
        bool isAtSurface = node is InletNode;
        if (node.DownstreamConveyance != null && !node.DownstreamConveyance.IsSubSurface)
          isAtSurface = true;
        foreach (Conveyance upConveyance in node.UpstreamConveyances)
          if (!upConveyance.IsSubSurface)
            isAtSurface = true;

        if (isAtSurface)
        {
          long index = -1;
          try
          {
            index = _GetRasterIndex(node.Point);
          }
          catch (IndexOutOfRangeException) { }

          if (index != -1)
          {            
            if (indexedSurfaceNodes.ContainsKey(index))
            {
              // In cases where multiple nodes occupy the same cell, transition nodes have priority
              if (node is TransitionNode)
                indexedSurfaceNodes[index] = node;

              string message = string.Format("Two surface nodes found in same raster pixel at ({0}, {1}). Prioritizing transition node if it exists.", node.Point.X, node.Point.Y);
              OnWarning(new WarningEventArgs(message));
            }
            else
            {
              indexedSurfaceNodes[index] = node;
            }
          }
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, "Nodes indexed."), true);
      return indexedSurfaceNodes;
    }

    private Dictionary<NodePoint, NodePoint> _SnapEndpoints(StormDrainNetwork overlandNetwork, Dictionary<long, NodePoint> indexedSurfaceNodes)
    {
      //TODO: Simplify any modified conveyance lines that are snapped here (to remove possible duplicate vertices)

      OnProgress(new ProgressEventArgs(0, 1, "Snapping surface flow lines to pipes..."), true);

      int progress = 0;
      int total = overlandNetwork.Nodes.Count;
      int increment = _CalcIncrement(total);

      Dictionary<NodePoint, NodePoint> nodeMapping = new Dictionary<NodePoint, NodePoint>();

      // We need to clip stream ends by a small amount because the StreamToFeature tool extends lines to the edge of the cell.
      double cellSize = ((IRasterAnalysisProps)_flowDirRaster).PixelWidth;
      double clipDistance = 0.01 * cellSize;
      IPoint indexPoint = new PointClass();

      //TODO: Track association b/w overland end point nodes and transition nodes from storm network
      foreach (NodePoint node in overlandNetwork.Nodes)
      {
        if (node.DownstreamConveyance == null && node.UpstreamConveyances.Count > 0)
        {
          // At downstream end
          IPolyline line = node.UpstreamConveyances.First().Line;
          line.QueryPoint(esriSegmentExtension.esriNoExtension, line.Length - clipDistance, false, indexPoint);

          long index = _GetRasterIndex(indexPoint);
          if (indexedSurfaceNodes.ContainsKey(index))
          {
            node.Point = indexedSurfaceNodes[index].Point;
            nodeMapping.Add(node, indexedSurfaceNodes[index]);
          }
          else
          {
            // No inlet or other node to snap to, but we want to pull point back from edge of cell.
            // Use clone because indexPoint is recycled.
            node.Point = (IPoint)((IClone)indexPoint).Clone();
          }

          // Match upstream conveyance lines to new node location
          foreach (Conveyance conveyance in node.UpstreamConveyances)
          {
            conveyance.Line.ToPoint = node.Point;
          }
        }
        else if (node.UpstreamConveyances.Count == 0 && node.DownstreamConveyance != null)
        {
          // At upstream end
          IPolyline line = node.DownstreamConveyance.Line;
          line.QueryPoint(esriSegmentExtension.esriNoExtension, clipDistance, false, indexPoint);

          long index = _GetRasterIndex(indexPoint);
          if (indexedSurfaceNodes.ContainsKey(index) && indexedSurfaceNodes[index] is OutletNode)
          {
            node.Point = indexedSurfaceNodes[index].Point;
            nodeMapping.Add(node, indexedSurfaceNodes[index]);
          }
          else
          {
            // No outlet to snap to, but we want to pull point back from edge of cell.
            // Use clone because indexPoint is recycled.
            node.Point = (IPoint)((IClone)indexPoint).Clone();
          }

          // Match downstream conveyance line to new node location
          node.DownstreamConveyance.Line.FromPoint = node.Point;
        }
        else
        {
          // middle of stream (can happen at outlets)
          IPolyline line = node.DownstreamConveyance.Line;

          long index = _GetRasterIndex(line.FromPoint);
          if (indexedSurfaceNodes.ContainsKey(index) && indexedSurfaceNodes[index] is OutletNode)
          {
            node.Point = indexedSurfaceNodes[index].Point;
            nodeMapping.Add(node, indexedSurfaceNodes[index]);
            node.DownstreamConveyance.Line.FromPoint = node.Point;
            foreach (Conveyance conveyance in node.UpstreamConveyances)
            {
              conveyance.Line.ToPoint = node.Point;
            }
          }
        }

        if (++progress % increment == 0)
        {
          OnProgress(new ProgressEventArgs(progress, total, "Snapping surface flow lines to pipes..."));
        }
      }

      OnProgress(new ProgressEventArgs(1, 1, "Surface flow lines snapped."), true);

      return nodeMapping;
    }

    private void _ReplaceOverlandNetwork(StormDrainNetwork overlandNetwork, Dictionary<NodePoint, NodePoint> nodeMapping)
    {
      // Track existing overland flowlines to correct low flow lines that are absent from the new overlandNetwork
      ElementCollection<Conveyance> existingFlowlines = new ElementCollection<Conveyance>();
      foreach (FlowLine flowline in _stormNetwork.FlowLines)
      {
          existingFlowlines.Add(flowline);
      }

      OnProgress(new ProgressEventArgs(0, 1, "Adding new surface flow lines...", true), true);

      int progress = 0;
      int total = overlandNetwork.Conveyances.Count;
      int increment = _CalcIncrement(total);
      bool cancel = false;

      List<Conveyance> tempList = overlandNetwork.Conveyances.ToList();

      // Add new overland flow lines
      for (int i = overlandNetwork.Conveyances.Count - 1; i >= 0; i--)
      {
        Conveyance conveyance = tempList[i];

        //TODO: Consider changing to factory pattern for nodes/conveyances of a network to avoid id conflicts

        // We need to clone conveyance because of id conflicts between networks
        FlowLine flowline = new FlowLine();
        flowline.Line = conveyance.Line;
        flowline.DataSource = conveyance.DataSource;

        if (nodeMapping.ContainsKey(conveyance.DownstreamNode))
        {
          ////// Remove existing upstream flowlines (removes terminal flow lines)
          ////if (nodeMapping[conveyance.DownstreamNode].ID < 0)
          ////{
          ////  foreach (Conveyance existingFlowline in nodeMapping[conveyance.DownstreamNode].UpstreamConveyances.ToArray())
          ////  {
          ////    if (terminalLines.Contains(existingFlowline))
          ////    {
          ////      terminalLines.Remove(existingFlowline);
          ////      _stormNetwork.Conveyances.Remove(existingFlowline);
          ////    }
          ////  }
          ////}

          flowline.DownstreamNode = nodeMapping[conveyance.DownstreamNode];
        }

        if (nodeMapping.ContainsKey(conveyance.UpstreamNode))
        {
          // Remove existing flowline if it exists (removes terminal flow lines)
          FlowLine existingConveyance = nodeMapping[conveyance.UpstreamNode].DownstreamConveyance as FlowLine;
          if (existingFlowlines != null)
          {
            _stormNetwork.Conveyances.Remove(existingConveyance);
            existingFlowlines.Remove(existingConveyance);
          }
          flowline.UpstreamNode = nodeMapping[conveyance.UpstreamNode];
        }

        // Adding flowline to storm network will create missing nodes.
        _stormNetwork.Conveyances.Add(flowline);

        flowline.UpstreamNode.TotalAccumulatedFlow = conveyance.MinAccumulatedFlow;

        // Remove references to free memory for garbage collection
        tempList.RemoveAt(i);
        overlandNetwork.Conveyances.Remove(conveyance);

        if (++progress % increment == 0)
        {
          ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Adding new surface flow lines...", true);
          OnProgress(eventArgs);
          if (eventArgs.Cancel)
          {
            cancel = true;
            break;
          }
        }
      }

      // Now we can snap the existing flowlines that were created in _LinkNodesOverland
      // since we have completed raster processing
      foreach (Conveyance existingFlowline in existingFlowlines)
      {
        existingFlowline.Line.FromPoint = existingFlowline.UpstreamNode.Point;
        existingFlowline.Line.ToPoint = existingFlowline.DownstreamNode.Point;
      }

      OnProgress(cancel
        ? new ProgressEventArgs(0, 0, "Rebuilding surface network canceled.")
        : new ProgressEventArgs(1, 1, "Adding new surface flow lines..."), true);
    }

    private void _LoadStormNetworkNodes(StormDrainNetwork stormNetwork)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Loading saved storm network nodes...", true), true);

      int progress = 0;
      int total = _nodesClass.FeatureCount(null);
      int increment = _CalcIncrement(total);

      int idField = _nodesClass.FindField(ID_FIELD);
      int typeField = _nodesClass.FindField(TYPE_FIELD);
      int enabledField = _nodesClass.FindField(ENABLED_FIELD);
      int sourceClassField = _nodesClass.FindField(SOURCECLASS_FIELD);
      int sourceOidField = _nodesClass.FindField(SOURCEOID_FIELD);
      int globalIdField = _nodesClass.FindField(GLOBALID_FIELD);
      int accumField = string.IsNullOrWhiteSpace(_flowAccumField) ? -1 : _nodesClass.FindField(_flowAccumField);
      int loadField = string.IsNullOrWhiteSpace(_loadAccumField) ? -1 : _nodesClass.FindField(_loadAccumField);
      int tagField = string.IsNullOrWhiteSpace(_tagField) ? -1 : _nodesClass.FindField(_tagField);

      IFeatureCursor cursor = null;
      try
      {
        cursor = _nodesClass.Search(null, false);
        IFeature nodeFeature = cursor.NextFeature();
        while (nodeFeature != null)
        {
          object dbValue;
          dbValue = nodeFeature.get_Value(idField);
          if (dbValue != DBNull.Value)
          {
            int id = (int)dbValue;

            dbValue = nodeFeature.get_Value(typeField);
            NodeType type = dbValue == DBNull.Value ? NodeType.Node : (NodeType)(int)dbValue;

            NodePoint node;
            if (type == NodeType.Inlet)
            {
              node = new InletNode(id);
            }
            else if (type == NodeType.Outlet)
            {
              node = new OutletNode(id);
            }
            else
            {
              node = new NodePoint(id);
            }

            node.Point = nodeFeature.ShapeCopy as IPoint;

            dbValue = nodeFeature.get_Value(enabledField);
            short enabled = dbValue == DBNull.Value ? (short)1 : Convert.ToInt16(dbValue);
            node.Enabled = enabled != (short)0;

            if (sourceClassField != -1 && sourceOidField != -1)
            {
              string sourceClass = nodeFeature.get_Value(sourceClassField) as string;

              dbValue = nodeFeature.get_Value(sourceOidField);
              int sourceOid = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);

              int globalId = -1;
              if (globalIdField != -1)
              {
                dbValue = nodeFeature.get_Value(globalIdField);
                globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
              }

              if (!string.IsNullOrEmpty(sourceClass) && sourceOid != -1)
              {
                node.DataSource = new FeatureIdentifier() { FeatureClassName = sourceClass, ObjectId = sourceOid , GlobalId = globalId};
              }
            }

            if (accumField != -1)
            {
              dbValue = nodeFeature.get_Value(accumField);
              node.TotalAccumulatedFlow = dbValue == DBNull.Value ? 0.0 : Convert.ToDouble(dbValue);
            }

            if (loadField != -1)
            {
              dbValue = nodeFeature.get_Value(loadField);
              node.TotalAccumulatedLoad = dbValue == DBNull.Value ? 0.0 : Convert.ToDouble(dbValue);
            }

            if (tagField != -1)
              node.Tag = nodeFeature.get_Value(tagField);

            stormNetwork.Nodes.Add(node);
          }

          if (++progress % increment == 0)
          {
            ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Loading saved storm network nodes...", true);
            OnProgress(eventArgs);
            if (eventArgs.Cancel)
              throw new UserCanceledException();
          }

          nodeFeature = cursor.NextFeature();
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} nodes loaded.", stormNetwork.Nodes.Count)), true);
    }

    private void _LoadStormNetworkConveyance(StormDrainNetwork stormNetwork)
    {
      OnProgress(new ProgressEventArgs(0, 1, "Loading saved storm network conveyance...", true), true);

      int progress = 0;
      int total = _conveyanceClass.FeatureCount(null);
      int increment = _CalcIncrement(total);

      int idField = _conveyanceClass.FindField(ID_FIELD);
      int upNodeField = _conveyanceClass.FindField(UPNODE_FIELD);
      int downNodeField = _conveyanceClass.FindField(DOWNNODE_FIELD);
      int enabledField = _conveyanceClass.FindField(ENABLED_FIELD);
      int subSurfaceField = _conveyanceClass.FindField(SUBSURFACE_FIELD);
      int sourceClassField = _conveyanceClass.FindField(SOURCECLASS_FIELD);
      int sourceOidField = _conveyanceClass.FindField(SOURCEOID_FIELD);
      int globalIdField = _conveyanceClass.FindField(GLOBALID_FIELD);
      int tagField = string.IsNullOrEmpty(_tagField) ? -1 : _conveyanceClass.FindField(_tagField);

      IFeatureCursor cursor = null;
      try
      {
        cursor = _conveyanceClass.Search(null, false);
        IFeature conveyanceFeature = cursor.NextFeature();
        while (conveyanceFeature != null)
        {
          object dbValue;
          dbValue = conveyanceFeature.get_Value(idField);
          if (dbValue != DBNull.Value)
          {
            int id = (int)dbValue;

            dbValue = conveyanceFeature.get_Value(subSurfaceField);
            short subSurface = dbValue == DBNull.Value ? (short)1 : Convert.ToInt16(dbValue);

            Conveyance conveyance;
            if (subSurface == (short)0)
            {
              conveyance = new FlowLine(id);
            }
            else
            {
              conveyance = new Conveyance(id);
              conveyance.IsSubSurface = true;
            }

            conveyance.Line = conveyanceFeature.ShapeCopy as IPolyline;

            dbValue = conveyanceFeature.get_Value(enabledField);
            short enabled = dbValue == DBNull.Value ? (short)1 : Convert.ToInt16(dbValue);
            conveyance.Enabled = enabled != (short)0;

            dbValue = conveyanceFeature.get_Value(upNodeField);
            int upNodeId = dbValue == DBNull.Value ? 0 : Convert.ToInt32(dbValue);
            if (stormNetwork.Nodes.Contains(upNodeId))
            {
              conveyance.UpstreamNode = stormNetwork.Nodes[upNodeId];
            }

            dbValue = conveyanceFeature.get_Value(downNodeField);
            int downNodeId = dbValue == DBNull.Value ? 0 : Convert.ToInt32(dbValue);
            if (stormNetwork.Nodes.Contains(downNodeId))
            {
              conveyance.DownstreamNode = stormNetwork.Nodes[downNodeId];
            }

            if (sourceClassField != -1 && sourceOidField != -1)
            {
              string sourceClass = conveyanceFeature.get_Value(sourceClassField) as string;
              dbValue = conveyanceFeature.get_Value(sourceOidField);
              int sourceOid = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);

              int globalId = -1;
              if (globalIdField != -1)
              {
                dbValue = conveyanceFeature.get_Value(globalIdField);
                globalId = dbValue == DBNull.Value ? -1 : Convert.ToInt32(dbValue);
              }

              if (!string.IsNullOrEmpty(sourceClass) && sourceOid != -1)
              {
                conveyance.DataSource = new FeatureIdentifier() { FeatureClassName = sourceClass, ObjectId = sourceOid, GlobalId = globalId };
              }
            }

            if (tagField != -1)
              conveyance.Tag = conveyanceFeature.get_Value(tagField);

            stormNetwork.Conveyances.Add(conveyance);
          }

          if (++progress % increment == 0)
          {
            ProgressEventArgs eventArgs = new ProgressEventArgs(progress, total, "Loading saved storm network conveyance...", true);
            OnProgress(eventArgs);
            if (eventArgs.Cancel)
              throw new UserCanceledException();
          }

          conveyanceFeature = cursor.NextFeature();
        }
      }
      finally
      {
        ComUtility.ReleaseComObject(cursor);
      }

      OnProgress(new ProgressEventArgs(1, 1, string.Format("{0} conveyance elements loaded.", stormNetwork.Conveyances.Count)), true);
    }

    /// <summary>
    /// Sets the output environment for this raster tool, using the specified extent for output.
    /// </summary>
    /// <param name="pRasterTool">The tool to set the ouput environment for.</param>
    /// <param name="extentProvider">An IEnvelope or RasterDataset that defines the output extent.</param>
    /// <remarks>The extentProvider determines the extent, but the default raster (the input flow direction grid) determines cell size
    /// and snap registration of the cells.</remarks>
    internal void SetAnalysisEnvironment(IRasterAnalysisEnvironment rasterTool, object extentProvider)
    {
      object inputGrid = _flowDirRaster;

      rasterTool.DefaultOutputRasterPrefix = "StrmTlsTmp";
      rasterTool.OutWorkspace = GetTempRasterWorkspace();
      rasterTool.OutSpatialReference = ((IGeoDataset)inputGrid).SpatialReference;
      rasterTool.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref inputGrid);
      rasterTool.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref extentProvider, ref inputGrid);
    }

    private ISpatialReference _GetPreferedSpatialReference()
    {
      ISpatialReference spatialRef = null;
      if (_geometricNetwork != null)
      {
        spatialRef = ((IGeoDataset)_geometricNetwork.FeatureDataset).SpatialReference;
      }
      else if (_flowDirRaster != null)
      {
        spatialRef = ((IGeoDataset)_flowDirRaster).SpatialReference;
      }
      return spatialRef;
    }

    internal IWorkspace GetTempRasterWorkspace()
    {
      IWorkspaceFactory pWSFactory = new RasterWorkspaceFactoryClass();
      return pWSFactory.OpenFromFile(_GetNewTempDir(), 0);
    }

    private string _GetNewTempDir()
    {
      string documents = Environment.GetFolderPath(Environment.SpecialFolder.MyDocuments);
      string tempDir = System.IO.Path.Combine(documents, "temp");
      tempDir = System.IO.Path.Combine(tempDir, Guid.NewGuid().ToString("D"));
      Directory.CreateDirectory(tempDir);
      return tempDir;
    }

    private long _GetRasterIndex(IPoint mapPoint)
    {
      return _GetRasterIndex(mapPoint.X, mapPoint.Y);
    }

    private long _GetRasterIndex(double x, double y)
    {
      int col, row;
      ((IRaster2)_flowDirRaster).MapToPixel(x, y, out col, out row);

      long c = col;
      long r = row;
      long w = ((IRasterProps)_flowDirRaster).Width;
      long h = ((IRasterProps)_flowDirRaster).Height;

      if (c >= w || c < 0 || r >= h || r < 0)
        throw new IndexOutOfRangeException("The point lies outside of the raster.");

      return r * w + c;
    }

    private INetworkFeature _GetNetworkFeature(int eid, esriElementType elementType, INetElements netElements)
    {
      // Check enabled field in Network Feature since it might not be in synch with network element enabled property
      int classID, userID, subID;
      netElements.QueryIDs(eid, elementType, out classID, out userID, out subID);

      INetElementDescriptionEdit elementDescription = new NetElementDescriptionClass();
      elementDescription.ElementType_2 = elementType;
      elementDescription.UserClassID_2 = classID;
      elementDescription.UserID_2 = userID;
      elementDescription.UserSubID_2 = subID;

      return _geometricNetwork.get_NetworkFeature(elementDescription);
    }

    private bool _IsDisabled(int eid, esriElementType elementType, INetElements netElements)
    {
      INetworkFeature netFeature = _GetNetworkFeature(eid, elementType, netElements);
      return _IsDisabled(netFeature);
    }

    private bool _IsDisabled(INetworkFeature netFeature)
    {
      if (netFeature == null)
        return true;

      bool disabled = false;
      IFeature feature = (IFeature)netFeature;
      int enabledField = feature.Fields.FindField("Enabled");
      if (enabledField > -1)
      {
        object enabledValue = feature.get_Value(enabledField);
        if (enabledValue != DBNull.Value)
        {
          disabled = !Convert.ToBoolean(enabledValue);
        }
      }
      return disabled;
    }

    private double _SafeConvertToDouble(object dbValue, double defaultValue = double.NaN)
    {
      if (dbValue != DBNull.Value)
      {
        try
        {
          return Convert.ToDouble(dbValue);
        }
        catch
        {
          return double.NaN;
        }
      }
      return double.NaN;
    }

    private int _CalcIncrement(int total)
    {
      int increment = total / 1000;
      if (increment == 0)
        increment = 1;
      return increment;
    }

    protected virtual void OnProgress(ProgressEventArgs eventArgs, bool log = false)
    {
      if (Progress != null)
      {
        Progress(this, eventArgs);
      }

      if (log)
      {
        Logger.Log(eventArgs.Message);
      }
    }

    protected virtual void OnWarning(WarningEventArgs eventArgs, bool log = true)
    {
      if (Warning != null)
      {
        Warning(this, eventArgs);
      }

      if (log)
      {
        if (eventArgs.InnerException == null)
        {
          Logger.Log(eventArgs.Message);
        }
        else
        {
          Logger.Log(string.Format("{0} ({1})", eventArgs.Message, eventArgs.InnerException.Message));
        }
      }
    }
  }
}
